#ifndef DYN_ODE_SIFS4_H
#define DYN_ODE_SIFS4_H

#include "dyn_ode_si.h"

/*! \brief Symplectic integrator with fixed stepsize of order 4.
 *
 * We implement the optimal fourth-order
 * symplectic integrator. We have taken
 * our values from Poulin (1997), which in turn
 * has taken them from McLachlan and Atela (1992).
 *
 * This integrator only works with Hamiltonians of the
 * form
 *  \f[ H(\mathbf{q},\mathbf{p}) = T(\mathbf{p})+V(\mathbf{q}), \f]
 * i.e., separable Hamiltonians.
 */

class SIFS4 : public SI
{
public:
    /*! Same as base constructor. */
    SIFS4(HamiltonianSystem& func,
          colvec _initCond,
          double _end,
          double _initStepsize)
        : SI(func,_initCond,_end,_initStepsize), M(4), order(4) {}

protected:
    /*! @name Coefficients of the expansion. */
    //@{
    const double c1 = 5.153528374311229364e-1;
    const double c2 = -8.5782019412973646e-2;
    const double c3 = 4.415830236164665242e-1;
    const double c4 = 1.288461583653841854e-1;
    const double d1 = 1.344961992774310892e-1;
    const double d2 = -2.248198030794208058e-1;
    const double d3 = 7.563200005156682911e-1;
    const double d4 = 3.340036032863214255e-1;

    const colvec c = {c1,c2,c3,c4};
    const colvec d = {d1,d2,d3,d4};
    //@}
};

#endif // DYN_ODE_SIFS4_H
